%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Explanation of what the file does
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% This function solves for the entire era (long timespan of T years)
% equilibrium of hats of price, actual labors, and wages, and levels of
% revenues and trade shares given labor supply and initial conditions in
% the counterfactual RUV economy
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Inputs needed
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Deficit matrix (DefEra), size I*T
% 2. Initial income by sector (labincsec), size I*S
% 3. Initial revenues (revbysec), size I*S
% 4. Initial trade shares (tradesharedat), size I*I*S
% 5. Hat trade costs (TAUmatC), size I*I*S*T
% 6. Hat technologies (TecHat), size I*S*T
% 7. Labor shares (labshares), size I*S
% 8. IO matrix (inoutmat), size I*S*S
% 9. Trade elasticities (Smat), size 1*S
% 10. Final good consumption shares (Amat), size I*S
% 11. Labor supply matrix (EllEra), size I*(S+1)*T
% 12. Delta DNWR matrix (DeltaMat), size I*S
% 13. Trade share levels in the baseline economy (ShaBa), size I*I*S*T
% 14. W dot in the baseline economy (WDBas), size I*S*(T-1)
% 15. L dot in the baseline economy (LDBas), size I*S*(T-1)
% 16. Parameter for updating AL algorithm (lambda), scalar
% 17. Tolerance parameter (tol), scalar
% 18. Agglomeration w.r.t. regional employment (AECR), flexible row vector
% 19. Agglome w.r.t. region-sector employment (AECRS), flexible row vector
% 20. Agglomeration parameter (aeatw), scalar, 1-> everyone, 2-> just USA
% 21. Place where non-USA regions start (M), scalar
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Outputs produced
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. Wages hat (Wmat), size I*S*(T-1)
% 2. Labor hat (Lmat), size I*S*(T-1)
% 3. Prices hat (PHma), size I*S*(T-1)
% 4. New revenue levels (RELE), size I*S*T
% 5. New trade share levels (TSLE), size I*I*S*T
% 6. Income by sector (Ytm1), size I*S*T
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%% Matlab functions invoked
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% 1. NECO_PeriodTatoSolver

function [Wmat,Lmat,PHma,RELE,TSLE,Ytm1]= NECO_Block2(DefEra,labincsec,...
revbysec,tradesharedat,TAUmatC,TecHat,labshares,inoutmat,Smat,...
Amat,EllEra,DeltaMat,ShaBa,WDBas,LDBas,lambda,tol,AECR,AECRS,aeatw,M)

% Define some useful variables

EllRel              = EllEra(:,2:end,:);
I                   = size(EllRel,1);
S                   = size(EllRel,2);
T                   = size(EllRel,3);
SDB                 = ShaBa(:,:,:,2:end)./ShaBa(:,:,:,1:end-1);
SDB(isnan(SDB))     = 0;

% Declare matrices of results

Wmat                = zeros(I,S,T-1);
Lmat                = zeros(I,S,T-1);
PHma                = zeros(I,S,T-1);

% Declare matrices of past levels

Ytm1                = zeros(I,S,T);
RELE                = zeros(I,S,T);
TSLE                = zeros(I,I,S,T);

% Fill some matrices in levels for the initial period (from the data)

Ytm1(:,:,1)         = labincsec;
RELE(:,:,1)         = revbysec;
TSLE(:,:,:,1)       = tradesharedat;

% Now run the time loop

for t=1:T-1

    Deficits        = DefEra(:,t+1);
    PastYnom        = Ytm1(:,:,t);
    PastSha         = TSLE(:,:,:,t).*SDB(:,:,:,t);
    
    if t==1
        PastCumLhat       = ones(I,S);
        TECmatrix         = TecHat(:,:,t);
    else
        PastCumLhat       = prod(Lmat(:,:,1:t-1).*LDBas(:,:,1:t-1),3);

        CumLhatB(:,:,1)   = ones(I,S);
        CumLhatB(:,:,2:t) = cumprod(LDBas(:,:,1:t-1),3);
        LLevB             = CumLhatB.*repmat(EllRel(:,:,1),1,1,t);
        LLevRegionB       = sum(LLevB,2);
        LRegDotB          = LLevRegionB(:,:,2:end)./LLevRegionB(:,:,1:end-1);
        
        CumLhat(:,:,1)    = ones(I,S);
        CumLhat(:,:,2:t)  = cumprod(Lmat(:,:,1:t-1).*LDBas(:,:,1:t-1),3);
        LLev              = CumLhat.*repmat(EllRel(:,:,1),1,1,t);
        LLevRegion        = sum(LLev,2);
        LRegDot           = LLevRegion(:,:,2:end)./LLevRegion(:,:,1:end-1);
        
        LRegHatrep        = repmat(LRegDot./LRegDotB,1,S,1); 
        sizem             = min(t-1,size(AECR,2));
        Relexponent       = repmat(permute(AECR( end-sizem+1:end),[1 3 2]),I,S,1);
        RaLRegDotrep      = LRegHatrep(:,:,end-sizem+1:end).^Relexponent;
        RelLmat           = Lmat(:,:,1:t-1);
        Relexponent2      = repmat(permute(AECRS(end-sizem+1:end),[1 3 2]),I,S,1);
        RaLDot            = RelLmat(:,:,end-sizem+1:end).^Relexponent2;
        TECmatrix         = TecHat(:,:,t) .* prod(RaLRegDotrep,3) .* prod(RaLDot,3);
    end
    
    if aeatw==2
        TECmatrix(M+1:end,:) = TecHat(M+1:end,:,t);
    end
    
    TAUmatrix       = TAUmatC(:,:,:,t);

    % Find the equilibrium for period t using a specially designed function

    Llim            = (EllRel(:,:,t+1)./EllRel(:,:,1)) ./ (PastCumLhat...
                      .*LDBas(:,:,t));
    Wlim            = DeltaMat./WDBas(:,:,t);
    maxiter         = 15000;
    [EqW,EqL,EqP,EqS,EqR,~] = NECO_PeriodTatoSolver(PastSha,TAUmatrix,...
    TECmatrix,labshares,inoutmat,Smat,Amat,PastYnom,Llim,Wlim,...
    Deficits,maxiter,tol,lambda,WDBas(:,:,t),LDBas(:,:,t));

    % Update some matrices like wages, employment, prices, revenues, etc

    Wmat(:,:,t)     = EqW;
    Lmat(:,:,t)     = EqL;
    PHma(:,:,t)     = EqP;
    Ytm1(:,:,t+1)   = EqW.*EqL.*PastYnom.*WDBas(:,:,t).*LDBas(:,:,t);
    RELE(:,:,t+1)   = EqR;
    TSLE(:,:,:,t+1) = EqS;

end

end
